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Abstract 

We discuss the confining and chiral-symmetry breaking properties of QCD with a large number 
of flavors Nf. In a Monte Carlo simulation of QCD with Nf = 16 staggered fermions, we find clear 
evidence of a first order bulk phase transition which separates phases with broken and unbroken 
chiral symmetry. This is consistent with extrapolations of earlier studies with smaller iVj, and is 
also as expected from general arguments. Thus, even when the perturbative renormalization group 
flow has a new infrared stable fixed point near the origin, lattice artifacts induce chiral symmetry 
breaking, and presumably confinement, at sufficiently strong coupling. 
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The confining properties of QCD are believed to depend sensitively on the number of light quarks 
species Nf. For example, when Nf > 17 the theory is not asymptotically free in the ultraviolet, and is 
presumably “trivial”, with free quarks and gluons. Chiral symmetry breaking is not expected to occur 
either. But an interesting change in the phase structure may occur even before Nf = 17. It was noted 
already in one of the first papers describing the two-loop beta function of QCD that the two-loop term 
changes sign at a value of Nf which differs from that of the one-loop term |]^ . This means that even for 
Nf < 17 a new infrared fixed point can occur in the renormalization group (RG) flow. From the two- 
loop beta function one obtains for SU{N), with Nf flavors of quarks in the fundamental representation 

ill 

{9*f ^ llN-2Nf 

167r2 -34N^+ l3NNf-3Nf/N 

for the fixed point. For SU{3) we thus get 


jg*)^ 16.5 - Nf 

167r2 -153 + l9Nf ’ 


where the right hand side of course should be positive. This leads to the constraint Nf > 8.053 (and 
Nf < 16.5). Taken at face value, this new fixed point thus occurs already for a number of flavors 
N* = 9. Of course, this prediction is perturbative, and hence can only be trusted if the new fixed point 
self-consistently occurs close to the origin, which is certainly not the case for Nf = 9. On the other 
hand, it is clear from these simple arguments that if we view the number of quark flavors Nf as a free, 
tunable (perhaps not even integer^) parameter, then the RG flow of QCD will be substantially altered 
as Nf is increased. This happens before asymptotic freedom is lost. However, it is clear that if Nf = 9, 
corresponding to ~ 0.417, we are so far from a perturbatively self-consistent region that the 

fixed point has no a priori meaning, and it is extremely dangerous to attempt to make much ado about 
Nf = 9. 

QCD with an infrared fixed point of the above kind was first analyzed in some detail by Banks and 
Zaks Ip, assuming that Nf was large enough to make the perturbative prediction valid. More recently, 
a study based on truncated Schwinger-Dyson equations has tried to determine more accurately the 
“critical number of flavors” N* at which the large-distance dynamics of QCD is governed by an infrared 
fixed point p. Another interesting approach, trying to go beyond the scheme-independent first two 
coefficients of the beta function and even beyond the (known) third coefficient in some given scheme, 
is based on the large-expansion [||. But it would obviously be desirable to consider the whole issue 
in a more stringent non-perturbative framework, in which the approximations can be controlled and 
improved in a systematic manner. This suggests the use of lattice gauge theory. 

We are therefore presently undertaking an investigation of the role played by Nf in lattice-regularized 
QCD, with Nf (continuum) species of staggered fermions, and the conventional Wilson gauge term 

(/3 = 6/ff2): 

Sgauge = f3Y.('^-\ReTvUp) (3) 

p 

where Up is the usual plaquette. 

In this letter we shall report on some results for Nf = 16. This is the largest integer value of Nf 
which is compatible with asymptotic freedom. On the other hand, it is also a value for which the 
prediction of a new perturbative infrared fixed point most likely can be trusted, since from eq. 
one gets {g *)‘^= 1/302. In principle, this fixed point could be singular, meaning that physical 
quantities cannot be obtained by just putting g = g*. However, in the following we shall assume 
that perturbation theory is asymptotically correct for small enough coupling as far as the first two 
(universal) terms in the /3-function are concerned. Then, for g* sufficiently small, there should be no 

^This is justified to any order of perturbation theory, since the Casimir operators are polynomials in Nf. 
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singular behaviour as the coupling g approaches g*. Close to the origin, the RG flow is thus restricted 
to lie in the interval 0 < g < g*, moving towards g* at large distances. From general arguments we do 
not expect this theory to be confining. What happens on the lattice? Consider the limit of infinitely 
strong bare gauge coupling; /3 = 0. For any number of quark flavors one can show that chiral symmetry 
is broken here i0- The effective Lagrangian, obtained by performing the gauge (link) integrations, 
can in this case explicitly be written in terms of meson and baryon fields Q. No free quark lines are 
permitted by the compact nature of the gauge integrations. It seems impossible not to associate this 
limit of extremely strong coupling with a confining phase. Since we have just argued from the other 
end of the coupling constant line that the theory is unconfining near the conventional continuum limit 
at g = 0, this means that a deconfining phase transition must occur somewhere in between. This phase 
transition is not expected to be associated with the infrared fixed point at g = g*. Rather, it should 
occur at much larger values of g, at the clash between the extreme strong-coupling lattice artifacts and 
the continuation beyond g* of the phase predicted by perturbation theory. 


We now present some of our results. We have performed simulations of the Nf = 16 theory in the 
range of gauge couplings 4.0 < (5 < 4.375. For our exploratory study of the phase structure we chose a 
relatively heavy staggered quark mass niq = 0.1. We then further explored the weak-coupling phase of 
the theory by varying the quark mass between 0.025 and 0.15. The R version of refreshed molecular 
dynamics algorithm was used for sampling of gauge configurations, with momenta refreshed every time 
unit. We used a microcanonical time step of 0.005, a conservative choice for which a convergence was 
observed in a multiflavor study by the Columbia group . 

We have measured, in addition to the chiral order parameter {'tp'ijj) and a variety of Wilson loops and 
correlation functions, the gauge-invariant fermion correlation functions appearing in the Fredenhagen- 
Marcu (FM) order parameter [10| (see below). We find strong evidence that the system undergoes a first 
order phase transition, with a metastability region ranging at least between /3 = 4.11 and /3 = 4.13. 
To clarify: our criterion for metastability is the existence of two phases at a given value of (3, with 
each phase having a lifetime of at least 200 microcanonical time units. In that region there are clear 
discontinuities in a number of observables: the chiral order parameter, the average action, and the string 
tension. Most importantly, as Figures |I]“|3| show, the location and the strength of the transition do not 
show any significant sensitivity to variations in the lattice size (although observation of metastability is 
difficult for small sizes). What we are seeing is thus most likely a bulk, rather than a finite-temperature 
phase transition. The nature of the transition is best revealed by the behavior of the chiral order 
parameter ('0V’) the string tension. The chiral order parameter decreases sharply in going from 
the low-/3 to the high-/? phase, suggesting that genuine chiral symmetry restoration occurs in the limit 
of vanishing quark mass (see below). Such restoration is usually also associated with a transition from 
confinement to deconfinement, but this is not a priori obvious. Of course, the notion of confinement 
in the presence of dynamical quarks in the fundamental representation of the gauge group is highly 
non-trivial (although an order parameter to distinguish between confinement and deconfinement in this 
case does exist, see below) due to the presence of quark-antiquark screening. Nevertheless, we would 
expect that the “string tension” formally defined, for intermediate distances, through a Creutz ratio of 
Wilson loops 

fw{x,t)w{x-i,t-iy 


a(x, t) = —In 


VVF(x,t-l)lR(a;-l,t) 


( 4 ) 


should drop discontinuously at a first-order deconfining phase transition. This is indeed what we observe 
in Figure The residual a on the relatively small coupling side (/? > 4.12) can at least qualitatively be 
understood as a “Coulomb”-like term. The contribution of a term of the form e/R (e is some effective 
coupling, here a free parameter) to a{x, t) is —e/x{x — l), so for x=3 we have a contribution —e/6 ~ 0.17 
if e is of order -1. The residual a on the relatively small coupling side (/? > 4.12) of Fig. can thus 
easily be understood by some effective “Coulomb” term, and not due to any string tension. A value 
e ~ — 1 is certainly not an unreasonable order of magnitude for a “Coulomb”-coefficient in the strong 
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coupling regime. For the much smaller number of flavors Nf = 2, a very accurate fit yields a value of 
e ~ —0.32 at /3 = 5.6, while for Nf = 0 (quenched simulations) the corresponding fit at the same value 
gives e ~ —0.29 [n.|. So the magnitude of this parameter increases with Nf at a fixed /3-value and is 
also expected to increase with the coupling g, consistent with our present estimate. Note also that the 
sign of this term is consistent with that of the Nf = 2 theory. 

A confinement criterion valid even in the presence of dynamical fermions has been proposed by 
Fredenhagen and Marcu [^]. We have tried to measure their corresponding order parameter 


P = 


lim 

R-^cc 


mO)U{R,T)i:iR)) 

^W{R,2T) 


R/T kept fixed 


( 5 ) 


as well. Here the numerator is a “staple” constructed by parallel displacement along a path of length 
T in the (arbitrarily chosen )“time” direction. The numerator is (the square root of) a Wilson loop 
of size R X 2T. The small sizes of our lattices obviously did not allow us to explore large space and 
time separations, as required by the FM order parameter. Moreover, our configuration sample was too 
small and did not permit us to determine the FM order parameter for the larger distances we actually 
had available. For smaller separations, where the measurement was possible, the hnite-distance version 
of order parameter (^) exhibited an somewhat unusual behavior. Namely, at zero separation, where 
the FM parameter coincides with the broken-phase value is larger than the symmetric-phase 

one. However, the FM parameter decreases with separation slower in the symmetric phase than it does 
in the broken phase, until finally the relative magnitude between the two is reversed. The only way 
to reconcile our data with the assumption, still unproven, that confinement implies chiral symmetry 
breaking, is by saying that the FM parameter approaches zero in the deconfined phase slower than it 
approaches a constant non-vanishing value in the confined phase. If so, considerably larger lattices are 
required for the FM parameter to be useful as a confinement criterion for these couplings, and with 
this number of Nf. 

At strong coupling, confinement/deconfinement can be checked explicitly by use of the FM order 
parameter. At infinitely strong coupling (/3 = 0) one can compute its value by either a hopping 
parameter expansion (in the case of Wilson fermions), or by an expansion in 1/m (in the case of 
staggered fermions). At infinitely strong coupling one finds, as expected, that the FM order parameter 
is nonvanishing (to leading order it is even constant, and the large-i? limit is reached immediately), 
indicating confinement. This holds for any number of flavors Nf. 


It is interesting to compare these results with earlier numerical work for smaller values of Nf, and, 
since the phase transition described here occurs at quite strong coupling, with existing strong coupling 
expansions. On the numerical side, zero-temperature chiral phase transitions for moderately large 
values of Nf were observed in refs. [Q, |^. For gauge group SU{3) and the same lattice action, they 
found a first-order phase transition at /3c = 4.47 ± 0.04 for Nf = 12 and lattice sizes 6“^ and 8^ [p^ . 
For the same gauge group, the first-order phase transition persisted for Nf = 8, and now occured at 
f3c — 4.85 on both 6^ and 8“^ lattices. More recently, in a very careful study, the Columbia group found 
a number of intriguing and puzzling features around the apparently bulk phase transition of the Nf = 8 
theory |^. Their value of the critical coupling equals /3c — 4.73 for lattices of sizes 16^ x 8 and 16^, at 
least in rough agreement with the above number. We note also the expected trend in the numbers: as 
Nf is increased, the phase transition moves toward stronger coupling. The fact that the phase transition 
is pushed into strong coupling as Nf grows, suggests that it eventually can be completely understood 
in terms of the strong coupling expansion. 

Although many of the gross features of lattice QCD with Nf = 16 thus seem to compare very well 
with earlier studies at smaller values of Nf, we nevertheless feel that a more precise characterization 
of the phase on the weak-coupling side of our observed transition is called for. As some steps in 
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Figure 1; The (3 dependence of the chiral order parameter for lattice sizes 8^ x 16 (diamonds), 12^ 
(triangles), and 6^ x 16 (squares). The quark mass is aniq = 0.1. The error bars are smaller than the 
plotting symbols. 



3.95 4 4.05 4.1 ^ 4.15 4.2 4.25 4.3 


Figure 2; The (3 dependence of the lattice action per site for lattice sizes 8^ x 16 (diamonds), 12^ 
(triangles), and 6^ x 16 (squares). The quark mass is aniq = 0.1. The error bars are smaller than the 
plotting symbols. 
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Figure 3: The f3 dependence of the string tension defined through a Creutz ratio 

— (hF(3,3)PF(2,2)) / (hF(3, 2)PF(2, 3)) for lattice sizes 8^ x 16 (diamonds), 12^ (triangles), and 6^ x 16 
(squares). 


this direction, we have investigated the behavior of the chiral condensate (^'0) and all the correlation 
functions as functions of the quark mass niq. The first question to answer is whether the weak-coupling 
side of the phase transition really is chirally symmetric in the limit of vanishing quark masses. To 
illustrate, we plot in Figure § as a function of aniq for the different /3-values quoted (at /3 = 4.125 
taking values from the weak-coupling phase). For all three /3-values the linear extrapolation of (^'0) 
from the smallest two quark masses to aniq = 0 is close to zero. (At /3 = 4.375 actually all three 
values lie on a straight line while, as discussed below, the aniq = 0.15, /3 = 4.25 is probably in the 
wrong phase.) This is our first evidence that the phase on the weak coupling side is chirally symmetric. 
Another piece of evidence comes from the fact that a quark-antiquark state with the pion quantum 
numbers in this phase does not behave as pseudo-Goldstone boson. This follows from Table 1, where 
we list a number of extracted mass values for the usual (lattice) particle spectrum. The pion mass 
definitely does not approach zero as niq is taken to zero, in this phase. It should be mentioned here 
that the fits required to extract the pion and rho masses on the chirally symmetric side were much 
more difficult to obtain, as compared with the chirally broken side (whereas the simplest exponential fit 
worked very well on the broken side, more terms were required on the symmetric). This could perhaps 
be an indication the whole bound-state particle concept may not be completely well-defined in the 
chirally symmetric phase. 

As niq increases, we expect that effectively fewer quark degrees of freedom participate in screening, 
and that hence the observed phase transition should move towards weaker coupling. Due to this, some 
(large) values of aniq on the presumed weak coupling side may actually correspond to the other phase. 
Consider figure where we plot the “string tension” a extracted from the Creutz ratios (^. For the 
point aniq = 0.15 at /3 = 4.25 we find a = 0.49(5) which is too large to appear in the figure. The point 
clearly belongs to the strong-coupling phase, as could already be suspected from the high value of {'ip'tp) 
in Figure Q. Otherwise the behavior of cr as a function of aniq is as expected: larger aniq corresponds 
to less screening, and hence, effectively, a larger value of a from (^). 
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Figure 4: The bare quark mass dependence of the chiral order parameter for the (3 values of 4.125 
(diamonds), 4.25 (squares), and 4.375 (triangles). The error bars are smaller than the plotting symbols. 
The dotted lines are drawn through the two lowest arug points at each value of (3. 


aniq 

/3 

TT 

P 

T^2 

P2 

N 

0.025 

4.25 

0.741 

0.910 

0.73 

0.91 

1.553 



±0.002 

±0.008 

±0.09 

±0.02 

±0.007 

0.05 

4.125 

0.749 

0.977 

0.92 

1.003 

1.605 



±0.004 

±0.005 

±0.01 

±0.007 

±0.007 

0.05 

4.25 

0.774 

0.972 

0.918 

1.012 

1.606 



±0.002 

±0.003 

±0.006 

±0.005 

±0.004 

0.05 

4.375 

0.780 

0.957 

0.892 

0.990 

1.586 



±0.025 

±0.004 

±0.008 

±0.003 

±0.005 

0.1 

4.0 

0.754 

1.57 

****** 

1.74 

2.37 



±0.001 

±0.03 

****** 

±0.025 

±0.02 

0.1 

4.125 

0.784 

1.47 

1.61 

1.64 

2.27 


broken 

±0.001 

±0.01 

±0.08 

±0.05 

±0.02 

0.1 

4.125 

0.845 

1.210 

1.181 

1.252 

1.90 


symmetric 

±0.002 

±0.008 

±0.006 

±0.006 

±0.01 

0.1 

4.25 

0.871 

1.123 

1.085 

1.152 

1.788 



±0.002 

±0.006 

±0.005 

±0.005 

±0.005 

0.1 

4.375 

0.882 

1.079 

1.039 

1.088 

1.737 



±0.002 

±0.004 

±0.005 

±0.007 

±0.003 

0.15 

4.25 

0.9323 

1.562 

1.676 

1.74 

2.51 



±0.0006 

±0.006 

±0.004 

±0.02 

±0.02 

0.15 

4.375 

0.992 

1.314 

1.289 

1.348 

1.989 



±0.002 

±0.014 

±0.006 

±0.005 

±0.0035 


Table 1: Summary of the hadron masses. 
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Figure 5: The bare quark mass dependence of the string tension for the /3 values of 4.125 (squares), 
4.25 (diamonds), and 4.375 (triangles). 


To understand these results in terms of the renormalization group, we have also made an attempt 
to extract a rough estimate of the beta function on the weak-coupling side of our phase transition. The 
rather crude method we have employed consists in measuring the ratio mp/m,r at one coupling (3 and 
for a given value of avriq. Choosing a different coupling (3' (on the same side of the phase transition), 
we vary aniq until the ratio mp/m,r matches the value at coupling f3. On an infinite volume, this gives 
us the (scheme dependent) scale change required to keep physics constant as we go from (3 to /?', and 
hence, indirectly, the beta function. These results indicate that the lattice spacing grows as we move 
towards weaker coupling, or, alternatively, that the beta function in the usual convention is positive in 
this region. This result is completely consistent with the two-loop perturbative beta function which 
precisely, for Nf = 16, passes through a zero close to the origin, and is positive beyond. If correct, 
the phase we are investigating on the weak coupling side is a pure lattice artifact, with no meaningful 
continuum limit (bounded as it is by a fixed point of vanishing correlation length on one side, and a 
first order phase transition on the other). The fact that this phase is strongly coupled but chirally 
symmetric (we call it strongly coupled since, like in ordinary QCD, the meson masses are much heavier 
than twice the quark mass) is from this point of view just another lattice artifact. 

If this picture were the whole story, we should presently be able to pursue directly our goal: to find 
the strict upper limit of light quark flavors Nf which is still compatible with chiral symmetry breaking 
and confinement. 0 Namely, we could lower Nf and follow the weak-coupling chirally symmetric phase. 
Measuring the beta function, we could find the precise point (as a function oi Nf) where the infrared 
fixed point disappears. However, there are a number of puzzling results, as well as a more fundamental 
problem, that prevent us from immediately doing so. For example, if we compare with the very precise 
study of the Nf = 8 theory in ref. [^, we note that both the phase transition and the weak-coupling 

^The weak-coupling phase we would see directly in such lattice simulations would correspond to the branch of the beta 
function beyond the infrared fixed point, with no continuum limit. But already the mere existence of an infrared fixed 
point at relatively weak conpling is presumable enough to prevent chiral symmetry breaking, and almost certainly strict 
confinement, also on the (physical) branch connected to the origin. 
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Figure 6: The bare quark mass dependence of the p to pion mass ratio for the P values of 4.125 
(diamonds), 4.25 (squares), and 4.375 (triangles). 


side of it look remarkably similar to what we find in the Nf = 16 theory. The apparently vanishing 
chiral condensate, the individual scaling of the mass spectrum in this phase: both are very similar in the 
Nf = 8 and the Nf = 16 theories. Naively, this would indicate that already a number Nf = 8 of quark 
flavors is incompatible with chiral symmetry breaking (and confinement). But in detail the results 
differ. For example, if we use the mass measurements in ref. Q to find the analogous beta function 
(by varying atriq, while keeping mp/rriT^ fixed, using mild extrapolations of their data), we find that it 
has precisely the opposite sign. Of course, the “beta function” we can extract in this manner is by no 
means universal, though taking the nucleon instead of the p gives qualitatively the same result. While 
it is a crude approximation to perform the matching in the same space-time volumes, the scale change 
between, for example P = 4.125, airiq = 0.05 and P = 4.25, airiq = 0.1 is only about 1.14 from the p 
mass, and hence the finite volume effects should not be too different and largely cancel. There seems 
to be a large anomalous scale factor for the quark mass. If the present “scheme” is sensible, we would 
not expect it to be able to introduce spurious zeros in the beta function so defined. Then the difference 
in sign of the beta functions for Nf = 8 and Nf = 16, which agrees with the perturbative prediction 
once, for Nf = 16, we are on the strong coupling side of the perturbative critical point, eq. (^), has 
to be taken seriously. We are presently investigating these scenarios in more detail. Another more 
fundamental problem concerns the non-universality of the beta function far away from fixed points 
of diverging correlation lengths. The precise location of an infrared fixed point at some intermediate 
coupling, and even the whole meaning of it, is scheme dependent. In lattice language, fixed points may 
develop which are to be viewed as no more than lattice artifacts. Therefore, even in principle it is a 
highly non-trivial task to pin-point the precise number of flavors at which confinement is lost. 


Finally we wish to mention that contrary to the arguments of ref. |14|, we see no reasons or 
numerical indications whatsoever for sensitivity to Nf on the extreme strong-coupling side. At very 
strong coupling there is no special significance to attach to the numbers Nf = 16 and Nf = 17, in 
contrast to the situation at weak coupling. Certainly, for Nf > 17 we do not expect to be able to 
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define a non-trivial continuum theory anywhere. But at strong coupling neither this value, nor values 
such as Nf = 6 or = 7 |14], play any special roles. The reason why the two pictures are totally 
compatible, is precisely to be found in the existence of a discontinuous bulk phase transition separating 
the two regimes. At least for staggered lattice fermions we thus find no support for the kind of scenario 
suggested, for Wilson fermions, in ref. [1T|. A more interesting question concerns the fate of the bulk 
phase transition as Nj is lowered, and, in particular, the nature of the phases on the weak-coupling 
side. Also this is presently under study. 
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